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A Short Tutorial on Inertial Navigation System and Global Positioning 

System Integration 


Kyle Smalling, Northrop Grumman 
Kenneth W. Eure, NASA Langley Research Center 


Abstract 

The purpose of this document is to describe a simple method of integrating Inertial Navigation 
System (INS) information with Global Positioning System (GPS) information for an improved estimate of 
vehicle attitude and position. A simple two dimensional (2D) case is considered. The attitude estimates 
are derived from sensor data and used in the estimation of vehicle position and velocity through dead 
reckoning within the INS. The INS estimates are updated with GPS estimates using a Kalman filter. This 
tutorial is intended for the novice user with a focus on bringing the reader from raw sensor measurements 
to an integrated position and attitude estimate. An application is given using a remotely controlled ground 
vehicle operating in assumed 2D environment. The theory is developed first followed by an illustrative 
example. 

1. Introduction 

With advancements in Micro-Electro-Mechanical Systems (MEMS) the estimate of vehicle 
position and attitude has become affordable and increasingly accurate as sensors have become smaller 
and cheaper. The integration of sensor outputs to achieve reasonable estimates remains an ongoing area 
of research. Integrating INS with GPS is most commonly done with tight coupling, which integrates sensor 
data at the measurement level, or loose coupling, which uses outputs of INS and GPS without knowledge 
of internal variables. The approach presented here is loose coupling. There is an abundance of literature 
pertaining to strapdown inertial navigation systems and integration with GPS [1], However, much of this 
literature is difficult to follow for the uninitiated. In order to greatly simplify complexity, the work 
presented here will be in two dimensions. Additional simplifications include a flat, nonrotating earth and 
short vehicle trajectories. These assumptions allow one to assume constant gravity, elevation, and 
magnetic field, all on a two-dimensional flat earth. 

The remainder of this paper is contained in the following 4 sections. Section 2 covers the 
mathematical development for attitude estimation around a single axis and INS-GPS fusion for derivation 
of a platform's 3 degrees-of-freedom (3DOF): (x nort/l ,x easf: ,iJ;). Next, Section 3 describes the 
experimental setup used to collect data for algorithm evaluation. Lastly, Sections 4 and 5 respectively 
present results and discuss conclusions. 

2. Mathematical development 

Reference frames play a large role in navigation. Before vehicle position and attitude can be 
estimated, one must know what the estimates are relative to. Here, we use two reference frames: the 
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vehicle body frame (x,y,z), and the navigation frame which is fixed to the earth (North, East, Down), (NED). 
In the body frame, we have the origin as the rotational center of the vehicle, x-axis out the front, and y- 
axis out the right side. The frame of the vehicle is assumed to always be coplanar with the north and east 
vectors; the z axis is normal to this plan pointing down. Since the distances traveled by the vehicle are 
very small in relation to the earth, we will assume a flat earth with perpendicular lines of latitude and 
longitude. Also, we will assume there are no variations in the local gravitational and magnetic fields. These 
assumptions together with loose coupling of the GPS-INS greatly simplify the mathematical development. 
Also, the development presented here is for a vehicle that is assumed to be traveling in 2D. This is a 
reasonable assumption since the experimental vehicle trajectory was over a localized flat surface. 

Estimation of vehicle position and attitude will be done in the navigation frame. It is assumed here 
that the senor frame and body frame are co-located and aligned. In order to be used for state estimation, 
the acceleration measurements onboard the vehicle must be converted from the measured body frame 
to the navigation frame. This may be accomplished by the simple matrix operations shown below in 
Algorithm 1 [2]. 

Algorithm 1 
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Figure 1 Diagram of two-dimensional vehicle body frame and navigation frame. Body frame is front of vehicle and right side (x, 

y) and navigation frame is North and East (N, E). 
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Figure 1 shows the vehicle body frame x and y and the navigation frame of north and east. The 
angular displacement is the vehicle yaw, ij;, and the angular rate is a), following the right hand rule with z 
pointing down. In this algorithm the outputs from the two axis magnetometer, m x (k) and m y (/c), are 
normalized by division of the magnetic vector magnitude, |m|. The yaw angle is found using the MATLAB 
atan2 function which preserves the correct quadrant. Care should be taken that a negative sign is placed 
on the magnetic y measure to match yaw as defined in the figure. The Directional Cosine Matrix (DCM), 
C{J, may then be formed using ij;. The DCM may then be used to compute the accelerations in the 
navigation frame. Throughout this algorithm, k denotes the integer time step, indicating that this 
algorithm is carried out every time step. The computed accelerations in the navigation frame, f n , may 
then be used to update the state vector in the Kalman filter, later described in Algorithm 3. Also, the yaw 
angle is combined with the gyro measures in the complementary filter to form an improved estimate of 
the vehicle attitude in Algorithm 2. 


There are numerous schemes to estimate vehicle attitude based on sensor data. One common 
technique is the implementation of a Bayesian estimator such as an Extended Kalman Filter (EKF), an 
Unscented Kalman Filter (UKF), or a particle filter (PF). Because of its simplicity, stability, and accuracy, a 
special form of the complementary filter developed by Mahony [3] will be used. This algorithm is 
computationally efficient and well suited for real time processing using limited computational resources; 
it is shown below in Algorithm 2. 


Algorithm 2 
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1 Normalization of the magnetic measurements is not necessary; however it may improve noise immunity. 
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C b (/c, k + 1) = /- 
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In Algorithm 2, the matrix <o b contains the output of the z-axis gyro sensor in the body frame. The 
matrix is that formed in Algorithm 1. The matrix C{J is an estimate of the true Directional Cosine Matrix 
formed by combining of Algorithm 1 with the gyro reading at time step k. T s is the sample rate. From 

CJJ the yaw is computed every time step. The subscripts on CJJ 2 and Cj^ represent the second row, first 

column and the first row, first column values respectively. The equation for C^(k, k + 1) forms the time 
step update. If a simpler form of integration is desired, one can use Euler integration with little loss of 
accuracy. A more accurate estimation of the acceleration in the navigation frame is obtained by using 
Algorithm 2 compared to Algorithm 1 since both magnetometer and gyro inputs are considered. The gain, 
K p , is set to one or greater and determines the relative contributions of the gyro and DCM to the attitude 
solution. In summary the inputs to Algorithm 2 are the gyro readings and DCM from Algorithm 1. The 
outputs are the attitude and an improved estimate of the DCM. The matrix CjJ is used to transform the 
accelerations, which are measured in the body frame, to the navigation frame. Once in the navigation 
frame, the accelerations may be used as inputs into Algorithm 3 to help estimate vehicle displacement 
and velocity. 

The acceleration expressed in the navigation frame may be combined with GPS measurements of 
position and velocity to perform state estimation of the vehicle in the navigation frame. This is done with 
a Kalman filter. In order to construct the Kalman filter, the vehicle dynamic equations are written in state 
space form. The state vector consist of the displacement, velocity, and acceleration in each component 
direction; that is north and east. 

x(k) = x(k - 1) + x(k - 1) ■ T s + ^ ■ x(k - 1) ■ T 2 
x(k) = x(k — 1) + x(k — 1) ■ T s 
x(k) = x(k — 1) | 2 

In matrix form we have 

x(k)l h T s l/ 2 ' T s 2 l [ x ( k_1 ) 

x ( k ) =01 T s ■ x(k-l) 

x ( k )J Lo 0 1 |x(k — 1) 

The above equation may be written in compact notation. 

x(k) = <D ■ x(k - 1) 

2 The acceleration is assumed not to change over one time step. 
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Here x is the state vector of displacement, velocity, and acceleration. 

In state space form, the Kalman filter may be expressed in orthogonal north and east directions as in 
Algorithm 3. Here, two channels are used to independently evaluate the state estimates, north and east. 

Algorithm 3 
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In Algorithm 3, the state to be estimated is x, and the measurements LatGps and LongGps are the latitude 
and longitude in meters displaced from an arbitrary starting point. These values are derived from the GPS 
readings 3 . Likewise VeIGps represents the respective velocities as measured from GPS. The matrices Q 
and R are the process and measurement noise covariance, respectively. In practice Q and R are tuned for 
desired performance. Algorithm 3 is an iterative process; it is executed every time a new GPS reading is 
acquired. It consist of two steps: a prediction step (indicated by a minus superscript), which computes an 
estimate of the state x before the latest GPS measurements are used, and an update step which uses the 
latest GPS and acceleration measurements. An error term, As, is formed between the measured and the 
predicted state. The Kalman gain, K, informs the algorithm how much weight to place on the error update. 
This weighted error term is used to update the state estimation. Since this process is repeated every time 
a new GPS reading arrives, the algorithm tracks the evolution of the state through time. 

3. Experimental Setup 

Algorithms 1 to 3 were used to post-process data obtained from a small radio controlled (R/C) car 
shown in Figure 2. The car was driven in loops around a local area and data from the INS and GPS were 
recorded. The purpose of this experiment was to demonstrate the use of the algorithms and explore 
implementation issues for practical vehicle attitude and position estimation. 


3 Because of the local flat earth approximation, the conversion from degrees of Latitude and Longitude is scalar. 
For our location in Hampton, VA we have 110,980 meters per degree Latitude and 88,892 meters per degree 
Longitude. 
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Figure 2. Radio controlled car used to gather data for navigation algorithms. Red circuit card is the VN200. 


The R/C car was an HPI Trophy Truggy Flux (part #:107018). A modification made to the vehicle 
was the replacement of the factory plastic body and wing with pieces of plywood in order to have an easy 
surface to mount the electronics. Suspension modifications were made to accommodate the extra weight 
from additional batteries and electronics. The factory 2-channel 2.4GHz radio was replaced with a 7- 
channel Spektrum DX-7s; this allowed additional channels to be used in conjunction with the autopilot on 
board, which was not used for this experiment. The factory off-road tires and rims were replaced with 
slightly smaller but much smoother and more rigid street tires in order to increase traction on the paved 
surfaces. With all batteries and modifications, the R/C car weighed 12.45 lbs., which is approximately 2 
lbs. heavier than stock. 


For this experiment a BeagleBoard Xm and a VN-200 GPS/INS were used to collect data. The VN- 
200 was configured to output a binary message stream which was saved to a USB flash drive for post- 
processing in MATLAB. To reduce this noise in the accelerometer data, a layer of Dubro 1/4 inch latex 
foam rubber was mounted in between the VN-200 and BeagleBoard. Additional noise reduction was 
realized by using a noise reduction filter within the VN-200 architecture. 
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4. Experimental Results 

The data recorded from the radio controlled car was post-processed to obtain estimates of vehicle 
attitude and position. Estimates were compared with those given by the VN200. Figure 3 shows the 
attitude of the car during a segment of the test. Since the test car was driven on a flat track, only yaw is 
shown. Pitch and roll are assumed to be negligible. 
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Figure 3. Yaw angle in degrees. Blue is from Directional Cosine Matrix, red is from VN200, and green is from complementary 

filter. 


7 


200 r 


Yaw 


0 


H* V 


DCM 
VN200 
Comp 


-200 


-400 


-600 


-800 


-1000 


0 .*1 

A 

s* V t 


. — 


-1200 


0 


50 


100 150 200 

Time in Seconds 


250 300 


350 


Figure 4. Yaw angle in degrees using unwrap 

As can be seen in Figure 3 and Figure 4, the estimates obtained from the DCM using the 
magnetometer only and those of the complementary filter closely follow the yaw computed by the VN200. 
Part of the difference may be attributed to the fact that the VN200 uses true north as the yaw reference 
while the estimates of the DCM and complementary filter use magnetic north. This would account for a 
difference of approximately 10°42' for the test track location. A zoomed in section is shown in Figure 5. 
Here, one can see by comparing the DCM solution with the complementary filter that the complementary 
filter is combining the magnetometer and gyro measures to produce a smooth attitude estimate. 
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Figure 5. Zoomed in section of yaw estimation showing smoothing action of complementary filter (green line). 

The attitude information is used by the inertial navigation system as described in Algorithm 3 to 
compute estimates of vehicle position. Vehicle position was computed in two independent channels using 
a Kalman filter to integrate the dynamic model with GPS readings to obtain a smooth position estimate. 
The plots shown in Figure 6 and Figure 7 are for the eastward velocity and position and those of Figure 8 
and Figure 9 are northward. 
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Figure 6. Eastward Velocity 



Figure 7 . Eastward Position 
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Figure 9. Northward Position. 
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Figures 6-9 show the velocity and position tracking of the Kalman filter over a large portion of the 
test path. As can be seen, close tracking is maintained throughout. In order to see an advantage gained 
by incorporating INS with GPS, a zoomed in section of the vehicle location is shown in Figure 10 and Figure 
11 . 


Eastward Position, Longitude 



Figure 10. Figure illustrating improvement of Eastward position estimates, red, by Kalman filter. 

Northward Position, Lattitude 



Figure 11. Figure illustrating improvement of Northward position estimates, red, by Kalman filter. 
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Figure 10 and Figure 11 illustrate the ability of the Kalman filter to use the dynamic model to fill 
in position measurements between GPS updates. 

5. Conclusions 

This tutorial describes a simple method for the integration of INS and GPS under a restrictive set 
of assumptions. The assumptions of a flat nonrotating earth with constant fields brought much 
simplification to the navigation problem. The purpose of this document is to offer a simple introduction 
into the complex field of strapdown inertial navigation systems and GPS/INS integration. The basic 
principles covered here may be extended to more complex navigation and filtering problems. 


Appendix A 

In this appendix the Kalman filter used in Algorithm 3 will be derived from Bayes' theorem [4], 
Consider the state space equation show below. 

*k ~ AXk—i + C{k-1 

y k = Hx k + r k 

Here x k is the n-dimension state at time step k, A is the state transition matrix, y k is the vector of sensor 
measurements and H describes the connection between the states and measurements. The vectors 
q k _ 1 and r k are independent zero mean Gaussian process and measurement noise vectors respectively. 
Since the system described above is linear and excited by Gaussian noise, the state and measurement 
vectors must also be Gaussian. The probability density function, pdf, for a multivariate Gaussian 
distribution is shown below where we use the notation N (random variable \ mean, covarience ). 

1 / 1 _ \ 

N(x\m,P)= 7 - exp I — ( x — m) T P 1 (x — m) 

(2tt) /2\P\ V 2 V 2 ) 

In the above |P| is the determinant of the covariance matrix P and m is the mean state vector. The Kalman 
filter has two steps: prediction and update. In the prediction an estimate of the state at time step k is 
computed based on past measurements y 1;fe _ 1 before measurement y k is taken. In the update the 
measurement y k is used to update the state estimation. The prediction and update steps occur at every 
time step making the Kalman filter a recursive algorithm. The recursive formulas are derived below. 

For the prediction at time step k, consider the joint pdf shown below conditioned on past 
measurements y^-i, here the equality follows from conditional probability. 

P(x k ,x k - 1 I yi:fc-l) = p{x k I * fc _ 1 ,y 1 :k _ 1 )p(x k _ 1 lyi:k-l) 

Because the first term on the right hand side has no dependence on the measurements y, we may write 
the above equation as shown below. 

P(x k ,x k _ i | y-i-k-i ) = p(x k | x fe _ 1 )p(x fe _ 1 | y 1± _ J 
Both sides may be integrated over x k _ x to obtain the Chapman-Kolmogorov equation. 
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PUfc I *k-i )pUk-i I yi:k-i ) daCk-1 


/> 00 n 00 

p(x fe ,x fe _i | yi; fe _i)dx fe _i = I 

J — CO J — 00 

The joint distribution now becomes the marginal distribution conditioned on past measurements yi ;k _i. 
One can solve this equation to obtain the desired pdf forthe predicted state mean mf and covariance Pjf . 

n 00 

PUfc I yi:fc-l ) = p( X k | X k _! )p( X k _! | y 1;fe _ 1 ) dx fe _ 1 

j — CO 

The pd/s in the above equation are Gaussian and may be expressed as shown below. 

x k I ™-k ’ p k ) = I W(x k I ^x k -i, Q )N( x k _! | m k _ 1 , P k _! )dx k _ 1 

J — CO 


Here Q is the process noise convenience matrix. Given the past values of m k _ 1 and P k _ 1 from the previous 
time step the above equation may be solved for the current predicted mean and covariance values 
m k and Pjf. These values are estimates based on measurements up to and including k-1, but do not 
include the measurement y k . The solution is shown below 4 . 

mf = Am k _! 

Pk = AP k ^A T + Q 

The above equations are predictions of the mean of the state x k and covariance at time step k, but without 
using the measurement y k . This is to say all measurements before k were used to give a best estimate of 
the state and covariance based on past measurements, state estimations and knowledge of the state and 
measurement distributions. The next step is to update this estimate using the information obtained in 
measurement y k . 

Bayes' rule is used in the update step. Recall Bayes' as shown below for random variables A and 
B. 


p(A IB) = 


p( B | A )p(A) 
p(B) 


In general Bayes' rule can also be expressed with an additional conditional variable C as shown below. 

p(B | A, C )p(A | C) 


p(A | B,C) = 


p(B \ C) 


For our derivation, we can use Bayes' rule as shown below. Here we have an estimation of x k that 
includes y k for the update step. 


pO/c I yk.yi-k-i ) = 


p(y/c I *fc,yi:fc-i)pQfc I yi:k-i) 
p(yfe I yi:k-i) 


4 For the Gaussian linear system case the Chapman-Kolmogorov equation may be explicitly solved. Alternatively 
note that the expected value E{x fc } = E{Ax k _ 1 + q k --i} = Am k _ 1 . Also the covariance cov{i4x fc _ 1 + q k ^-i} = 
E{(Ax k ^ + q k - Am k - 1 )(.Ax k - 1 + q k - x - Am^Y] = AE^x^xl^ - + Q = AP k _ 1 A T + Q 

where E is the expectation operator and the time dependence of Q has been suppressed. 
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Since our derivation involves Gaussian distributions, we are only interested in computing the mean and 
covariance of the state update. To this end we will only examine the exponent of the distribution; which 
contains the mean and covariance, x k and P k respectively. We can express this as the proportion shown 
below after recognizing that p{y k \ x k ,y 1:k _ 1 ) = p(y k I x k ) since there is no dependency on y\±-i in 
our state space model. Also we have written y k ,y\- k -\ in compact notation as y 1:k . 

POfc lyi:k) K P(7k \x k )p(x k I y 1:k -i) 

The pdfs in the right hand side of this equation are known. The first term is simply the measurement 
equation pdf and the second is that determined in the previous prediction step. We may now express the 
update step in terms of the respective Gaussian distributions. 

N{x k | m k ,P k ) oc N(y k \ Hx k ,R)N(x k \ m k ,P k ) 

Here R is the measurement noise covariance matrix and mf,P k are the estimates computed in the 
prediction step. The goal is to find m k ,P k which are the updated mean and covariance based on the 
current measurement y k . One can do this by completing the square in the product of the two Gaussian 
pdfs given above. 

N(x k | m k , P k ) oc 

- — .n/ 1 , ,1/ e xp[-\{y k -HXkYR- 1 ^ -Hx k )) — jr ex P \ ~\(x k - mfY^y 1 ^ - mf ) ) 

(2tt) /2|B|V 2 V 2 /(2tt) /2|P-|V 2 V 2 J 

Consider the exponent only. The expansion is shown below. In this expansion, care must be taken to maintain the 
order since this is a matrix equation and multiplication does not commute. 

-\{yk R ~ 1 yk-y k R ~ lH x k -xlH T R- 1 y k + xlH T R- 1 Hx k + xl{p k y 1 x k -xl{p k y 1 m-f 
~ (rn k ) T (P k y 1 x k + (mfc) 7 ’(P fe “) _1 mfc) 

We need to shape the above expression into the form (x^ — m k ) T P k 1 (x k — m k ) in order to identify the 
updated mean m k and covariance P k . In order to do this we gather like terms and drop any term that 
does not contain x k . The dropped terms can be ignored since they do not contribute to the expression for 
the Gaussian mean and covariance. Factoring out x k and recalling that the transpose of a covariance 
matrix is just the matrix itself we obtain the expression below. 

~^{ x k{H T R~ 1 H + (P-)- 1 )x k - x k (H T R~ 1 y k + - (y^R^H + {mff {Pfy')x k } 

In order to make the math more tractable while completing the square we define some variables. 

V = ^R^H + ( P~y 1 
b = H T R~ 1 y k + (P”) _1 mfc 
Using the above, we can write our expression as shown. 

x kVx k — x k b — b T x k 

In order to complete the square we can add a term to this expression as long as it does not contain x k . 
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X k Vx k — x kb — b Tx k + b T V 1 b 

Also, recall that I/I/ -1 = I can be inserted into any term to help complete the square. 

xlVx k - x k VV~ 1 b - b T VV~ 1 x k + P T I/ -1 I/I/ -1 6 = (x k - V~ 1 b) T V(x k - I/ -1 P) 

The fact that the above equality holds can be verified by expanding the right side to get the left. We now 
have the desired form and can identify the updated mean and convenience matrices. 

m k = V-H = (H T R~ 1 H + (1 p-)- 1 ) -1 (// 7 ’ J R- 1 y fc + (P; c y lm k) 

p k = i/ -1 = [h t r~ 1 h + (P, - )- 1 ) -1 

To make the computations more efficient, we use the matrix inverse lemma shown below where A,U,C 
and V are simply placeholders to demonstrate the lemma. 

(. A + ucvy 1 = A- 1 - ^ -1 f/(C -1 + I M -1 £/) -1 I/,4 -1 
Pk = ((P^r 1 + H-'-R-'H)- 1 = P k - P k H T (R + //P k -// 7 ’) -1 //P k - 
m k = {P k - P k H T (R + HP k H T )~ 1 HP k }(// T P -1 y k + (P - ) -1 m k ) 
m k = Pk H T R~ 1 y k + m~ k - P k H T (R + HP k H T )~ 1 HP k H T R~ 1 y k - P k H T (R + HP k H T )~ 1 Hm k 
We need to simplify the y k term. 

Pfc/Z^P -1 - P k H T {R + HP k H T Y 1 HP k H T R- 1 = 

P ^// 7 ^ -1 - P k H T {R + HP k H T y 1 {R + HP k H T )R~ 1 + P k H T {R + //P k // T ) -1 = 

P k H T (R + HP k H 7 )* 1 

Replacing the y k coefficient with the above simplified expression one may now write m k as shown below. 
The previously derived prediction equations and covariance update are also shown for completeness. 

Prediction step: 

m k = Am k _! 

Pfc = AP k . ,A T + Q 

Update step: 

m k = m~ k + P k H T {R + HP k H r y\y k - Hm k ) 

Pk = Pk ~ PkH T (R + RP k H r )^HP k 

The equations are iterated every time step as in Algorithm 3. 

Appendix B 

The purpose of this section is to provide MATLAB code to implement Algorithms 2 and 3. The code 
snippets are partial and are not intended to be functioning code, but are presented only as an aid to the 
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reader. Both blocks of code run in a loop that executes every time step. The code for Algorithm 2 is shown 
below. 


mnorm = norm([mx(k) my(k)]); %mx(k) my(k) from magnetometer 
mx(k) = mx(k) /mnorm; my(k) = my (k) /mnorm; 
psi = atan2 ( -my ( k) , mx ( k) ) ; 

DCM = [cos (psi) -sin (psi); 

sin (psi) cos (psi)]; 

YAW(k) = 180/pi*psi; %YAW from magnetometer 
gyro = [ 0 -wz (k) ; %wz(k) from z axis gyro 

wz ( k) 0 ] ; 

omaga = 0.5*(Rh'*DCM - DCM ' *Rh) ; 

OM = gyro + kp*omaga; % kp >= 1, gain constant 
nOM = norm(OM); %For exact solution 

A = eye (2 ) - OM. *sin (nOM*Ts) . /nOM + (0M A 2) .* (1-cos (nOM*Ts) ) /nOM A 2; 

Rh = Rh*A ' ; 

yaw(k) = 180/pi*atan2 (Rh (2, 1) ,Rh (1, 1) ) ; %yaw. Complementary Filter 

The above code returns Rh which is an improved estimate of the DCM. From Rh the yaw angle is 
calculated. 

The code for Algorithm 3 is shown below. This block is for the longitude channel. The inputs are the time, 
GPS longitude, GPS velocity, and eastward acceleration in the navigation frame. 

Ts = Time (k) -Time ( k— 1 ) ; %Time step 
A = [1 Ts 0.5*Ts A 2; 0 1 Ts; 0 0 1]; 
x ( : , k) = A*x(:,k-1); %Predict 
P = A* P*A ' + Q; 

err ( : , k) = [LongGpsE ( k) ; VelNedE(k); accel^nav (2 , k) ] - x ( : , k) ; 

S = P + R; 

K = P*(S\I); 

x ( : , k) = x ( : , k) + K*err ( : , k) ; %Update 
P = (I - K) *P; 

The estimated state x ( : , k) combines both the INS and GPS measurements for improved estimates of 
position, velocity, and acceleration. The latitude channel follows in like manner. 
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